#### Simulation AR(1) and AR(2) for testing the first order autocorrelation
#### In sample properties
#### Generate Table 1 and 2
#### The user can change vecrho for the autocorrelation values and vecrho1 for the alternatives used for the roots
#### alpha, the significance level, vecn, the vector of sample sizes, nsimu the number of simulations drawn are also global parameters fixed by the user
rm(list=ls())

### Library #####
library(MASS)

###Rho
rho0=0.5
vecrho1=c(0,0.2,0.5)


## Significance level
alpha=0.05

##Sample sizes
### We estimate the parameters on R values and test on P values.
listep=c(50,100,150,250,500)
lister=c(50,100,150,250,500)

maxp=max(listep)
nnr=length(lister)
nnp=length(listep)

### Number of replications
nsimu=10000


for (rho1 in vecrho1){
rho11=rho0+rho1
rho12=-rho0*rho1
### Output File
nom=paste("Table3plus with rho=",rho0," and rho1=",rho1,".out")



##### Function which tests the first autocorrelation with parameter uncertainty
##### We compare robust moment and correction strategy, see West  and McCracken (1998)
##### It gives back a vector of ones and zeros depending on the test decision


testfirst<-function(r,p,rho11,rho12,nns){
set.seed(nns*10)
buffer=100
#n:taille fichier
n=r+p+buffer
eps=rnorm(n)
denvar= ((1-rho12)^2-rho11^2)*(1+rho12)/(1-rho12)
y0=rnorm(1)/sqrt(denvar)
ym1=rnorm(1)/sqrt(denvar)

ytsim=filter(eps,c(rho11,rho12),method="recursive",sides=1,init=c(y0,ym1))
yt=ytsim[(buffer+1):n]

##### Rolling estimator ######
rolling<-function(h){
yt0=c(ytsim[buffer],yt)
ytm1=yt0[h:(r+h-1)]
det=mean(ytm1^2)-(mean(ytm1))^2
hatmu=(mean(ytm1^2)*mean(yt[h:(r+h-1)])-mean(ytm1)*mean(ytm1*yt[h:(r+h-1)]))/det
hatrho=(-mean(ytm1)*mean(yt[h:(r+h-1)])+mean(ytm1*yt[h:(r+h-1)]))/det
epsthat=yt[h:(r+h-1)]-hatrho*ytm1-hatmu # Estimated inovation
epst=yt[h:(r+h-1)]-rho0*ytm1  # true innovation
epshatnew=yt[r+h]-hatrho*yt[r+h-1]-hatmu
epsnew=yt[r+h]-rho0*yt[r+h-1]
cor=(epshatnew*epsthat[r])/(mean(epsthat^2))
mt=(epshatnew*epsthat[r]-(1-hatrho^2)*epshatnew*(yt[r+h-1]-mean(ytm1)))/(abs(hatrho)*mean(epsthat^2))# our moment on estimated residuals
mtalt=(epshatnew*epsthat[r]-(1-hatrho^2)*epshatnew*(yt[r+h-1]-hatmu/(1-hatrho)))/(abs(hatrho)*mean(epsthat^2))# plug the mean
cor0=(epsnew*epst[r])/(mean(epst^2))
mt0=(epsnew*epst[r]-(1-rho0^2)*epsnew*yt[r+h-1])/(rho0*mean(epst^2))# our moment on true innovation
return(cbind(hatrho,cor,mt,cor0,mtalt,epshatnew,epsthat[r],yt[r+h-1]))
}
##### Fixed estimator ######
fixed<-function(h){
yt0=c(ytsim[buffer],yt)
ytm1=yt0[1:r]
det=mean(ytm1^2)-(mean(ytm1))^2
hatmu=(mean(ytm1^2)*mean(yt[1:r])-mean(ytm1)*mean(ytm1*yt[1:r]))/det
hatrho=(-mean(ytm1)*mean(yt[1:r])+mean(ytm1*yt[1:r]))/det
epst=yt[r+h]-hatrho*yt[r+h-1]-hatmu
epstm1=yt[r+h-1]-hatrho*yt[r+h-2]-hatmu
epsthat=yt[1:r]-hatrho*ytm1-hatmu
cor=(epst*epstm1)/(mean(epsthat^2))
mt=(epst*epstm1-(1-hatrho^2)*epst*(yt[r+h-1]-hatmu/(1-hatrho)))/(abs(hatrho)*mean(epsthat^2))# our moment on estimated residuals
mtalt=(epst*epstm1-(1-hatrho^2)*epst*(yt[r+h-1]-mean(ytm1)))/(abs(hatrho)*mean(epsthat^2))# our moment on estimated residuals
return(cbind(cor,mt,epst,epstm1,yt[r+h-1],mtalt))
}

matrolling=t(sapply(seq(1,p,1),rolling))
matfixed=t(sapply(seq(1,p,1),fixed))

##### test rolling West-McCracken #####
dep=as.matrix(matrolling[,6])
expl=as.matrix(cbind(matrolling[,7],rep(1,p),matrolling[,8]))
matv=solve((t(expl)%*%(expl))/p)%*%((t(expl)%*%((matrolling[,7]^2)*expl))/p)%*%solve((t(expl)%*%(expl))/p)
matalph=solve(t(expl)%*%(expl))%*%t(expl)%*%dep
testwm_rol=abs(sqrt(p)*matalph[1]/sqrt(matv[1,1]))

#### test fixed West-McCracken #####
dep=as.matrix(matfixed[,3])
expl=as.matrix(cbind(matfixed[,4],rep(1,p),matfixed[,5]))
matv=solve((t(expl)%*%(expl))/p)%*%((t(expl)%*%((matfixed[,4]^2)*expl))/p)%*%solve((t(expl)%*%(expl))/p)
matalph=solve(t(expl)%*%(expl))%*%t(expl)%*%dep
testwm_fix=abs(sqrt(p)*matalph[1]/sqrt(matv[1,1]))


### First test with true value #####
### Ours ##### 
xi4=sqrt(p)*mean(matrolling[,3])
### Ours ##### 
xi5=sqrt(p)*mean(matrolling[,5])
### Test with the rolling estimator #####
### Ours ##### 
xi6=sqrt(p)*mean(matrolling[,3])/sqrt(var(matrolling[,3]))
### Test with the fixed estimator #####
### Ours ##### 
xi8=sqrt(p)*mean(matfixed[,2])
### Ours ##### 
xi9=sqrt(p)*mean(matfixed[,6])
### Test with the rolling estimator #####
### Ours ##### 
xi10=sqrt(p)*mean(matfixed[,2])/sqrt(var(matfixed[,2]))

return((abs(c(xi4,xi6,xi8,xi10,testwm_rol,testwm_fix,xi5,xi9))>qnorm(1-alpha/2))*1)
}




####### We do it in parallel ####################################
library(foreach);
library(iterators);
library(doParallel)
nworkers <- detectCores()
cl <- makePSOCKcluster(nworkers-1)
#clusterSetRNGStream(cl, c(1,2,3,4,5,6,7,))
registerDoParallel(cl)

matsortie4=matrix(NA,nnr,nnp);
matsortie5=matrix(NA,nnr,nnp);
matsortie6=matrix(NA,nnr,nnp);
matsortie8=matrix(NA,nnr,nnp);
matsortie9=matrix(NA,nnr,nnp);
matsortie10=matrix(NA,nnr,nnp);
matsortie11=matrix(NA,nnr,nnp);
matsortie12=matrix(NA,nnr,nnp);

for (sr in 1:nnr){
	for (sp in 1:nnp){
	matrej<- foreach(ns=1:nsimu,.packages = c("nloptr","MASS"), .combine=rbind) %dopar% {testfirst(lister[sr],listep[sp],rho11,rho12,ns)}
	matsortie4[sr,sp]=mean(matrej[,1])*100
	matsortie5[sr,sp]=mean(matrej[,7])*100
matsortie6[sr,sp]=mean(matrej[,2])*100
matsortie8[sr,sp]=mean(matrej[,3])*100
matsortie9[sr,sp]=mean(matrej[,8])*100
matsortie10[sr,sp]=mean(matrej[,4])*100
matsortie11[sr,sp]=mean(matrej[,5])*100
matsortie12[sr,sp]=mean(matrej[,6])*100

}
}
stopCluster(cl)

####################################################################
####################################################################
library(xtable)
outputfile=paste(paste(nom,sep=""),".tex")
write("Table 3 + appendix, first autocorrelation test", file=outputfile,append=FALSE)
write(paste("rho0=",rho0),file=outputfile, append= TRUE)
write(paste("rho1=",rho1),file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)


####################################################################

write("ROLLING SCHEME",file=outputfile, append= TRUE)
write("**************",file=outputfile, append= TRUE)
rownames(matsortie4)=lister;
colnames(matsortie4)=listep;
print(xtable(matsortie4,caption="Est, rolling, robust"),digits=2,file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
rownames(matsortie5)=lister;
colnames(matsortie5)=listep;
print(xtable(matsortie5,caption="Est,rollingalt, robust"),digits=2,file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
rownames(matsortie6)=lister;
colnames(matsortie6)=listep;
print(xtable(matsortie6,caption="Est,rolling, robust,varinsample"),digits=2,file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
rownames(matsortie11)=lister;
colnames(matsortie11)=listep;
print(xtable(matsortie11,caption="Test rolling WM"),digits=2,file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)

####################################################################

write("FIXED SCHEME",file=outputfile, append= TRUE)
write("**************",file=outputfile, append= TRUE)
rownames(matsortie9)=lister;
colnames(matsortie9)=listep;
print(xtable(matsortie9,caption="Est,fixed,robust"),digits=2,file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
rownames(matsortie8)=lister;
colnames(matsortie8)=listep;
print(xtable(matsortie8,caption="Est,fixed,robustalt"),digits=2,file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)

rownames(matsortie10)=lister;
colnames(matsortie10)=listep;
print(xtable(matsortie10,caption="Est,fixed,robust,varinsample"),digits=2,file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
write("",file=outputfile, append= TRUE)
rownames(matsortie12)=lister;
colnames(matsortie12)=listep;
print(xtable(matsortie12,caption="Test fixed WM"),digits=2,file=outputfile, append= TRUE)
####################################################################
####################################################################
}